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The nonlinear frequency shift is derived in a transparent asymptotic form for intense Langmuir 
waves in general collisionless plasma. The formula describes both fluid and kinetic effects simultane¬ 
ously. The fluid nonlinearity is expressed, for the first time, through the plasma dielectric function, 
and the kinetic nonlinearity accounts for both smooth distributions and trapped-particle beams. 

Various known limiting scalings are reproduced as special cases. The calculation avoids differential 
equations and can be extended straightforwardly to other nonlinear plasma waves. 
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I. INTRODUCTION 

It has long been known that the nonlinear frequency 
shifts 6uj of collisionless plasma waves can depend on the 
wave amplitude a in various ways, even at small a. The 
contributions to Suj that stem from fluid effects typically 
scale as oc [1-6], whereas kinetic effects can produce 
contributions scaling as oc ^/a [7-21] or even decreasing 
with a [22-28]. The interest to such nonlinear dispersion 
relations (NDRs) has been revived recently in connec¬ 
tion with laser-plasma interactions [29] and the evolution 
of energetic-particle modes in tokamaks [30]. However, 
each of the existing theories that offers an explicit for¬ 
mula for Soj describes just one type of nonlinearity. The 
possible effect of multiple coexisting nonlinearities that 
yield different scalings in a was explored in Ref. [31], but 
only heuristically and for an artificial plasma model. The 
more general existing theories, such as that of Bernstein- 
Greene-Kruskal (BGK) waves [32], rely on formal solu¬ 
tions of the Vlasov-Maxwell equations and are not always 
easy to apply in practice due to their complexity. But is 
it possible to derive a comprehensive and yet tractable 
asymptotic NDR for general collisionless plasma? 

The answer is yes, and here we show how to do it 
explicitly. We build on the theory that was developed 
recently in Refs. [33-37] and represents a fully nonlin¬ 
ear kinetic version of Whitham’s average-Lagrangian ap¬ 
proach [38, 39]. (We assume that, when present, res¬ 
onant particles are phase-mixed, so the corresponding 
waves are of the BGK type [40].) This allows deriving 
the NDR directly from the wave Lagrangian, which is 
known, without solving any differential equations. In 
fact, all wave properties can be traced to the properties 
of a single function characterizing individual particles, 
namely, the normalized action j of a particle as a func¬ 
tion of its normalized energy r. We focus on Langmuir 
waves in one-dimensional electron plasma, but extending 
the theory to general waves is straightforward to do. 

Unlike in Refs. [34-37], where a related theory was 
constructed under the sinusoidal-wave approximation, we 
now allow for a nonzero amplitude of the second har¬ 
monic (higher-order harmonics are assumed negligible) 


and find this amplitude self-consistently. Kinetic and 
fluid nonlinearities are hence treated on the same foot¬ 
ing, and two types of distributions are studied in detail: 
(i) First, we consider particle distributions, /o, that are 
smooth enough near the resonance. We assume that a 
is small in this case; then kinetic nonlinearities scale as 
half-integer powers of a, and fluid nonlinearities scale as 
integer powers of a. In particular, it is shown for the 
first time that the fluid nonlinearity can be expressed in 
terms of the plasma dielectric function and, contrary to 
the common presumption, can have either sign, (ii) Next, 
we consider the effect of abrupt trapped-particle beam 
distributions superimposed on smooth /q. The small pa¬ 
rameter in this case is not a but rather the appropriately 
normalized trapped-particle density. We show that such 
beams contribute additively to the NDR. We revisit the 
case of flat-top trapped beams [26, 27], as an illustra¬ 
tion, and demonstrate that our asymptotic theory yields 
predictions for both the wave frequency and the second- 
harmonic amplitude that are virtually indistinguishable 
from those given by the exact solution of the Vlasov- 
Poisson system. Our theory then can be considered as 
advantageous over such solutions. This is because, while 
offering a reasonable precision, it is more transparent and 
flexible, i.e., can be used also when exact analytical for¬ 
mulas do not exist. 

Nonlinearities with both smooth particle distribution 
(including fluid and kinetic effects) and distributions with 
abrupt changes near the resonant region (beam distribu¬ 
tion) are treated on the same footing. In particular, we 
show, for the first time, that the fluid nonlinearity can 
be expressed in terms of the plasma dielectric function 
and, contrary to the common presumption, can have ei¬ 
ther sign. Nonlinearities due to abrupt distributions of 
trapped beams are described separately and can be added 
to the fluid and kinetic nonlinearities. 

The paper is organized as follows. In Sec. II, we briefly 
review the underlying variational approach and introduce 
basic notation. In Sec. Ill, we describe the wave model 
and asymptotics of some auxiliary functions derived from 
j(r). In Sec. IV, we derive the general expression for 
6u! for smooth /o and present examples, including the 
cases of cold, waterbag, kappa, and Maxwellian /q. In 


2 


Sec. V, we discuss the effect of trapped-particle beams 
superimposed on a smooth /g. In Sec. VI, we summarize 
the main results of our work. Some auxiliary calculations 
are also presented in appendixes. 

II. BASIC CONCEPTS AND NOTATION 
A. Wave Lagrangian density 

As reviewed in Refs. [33], the dynamics of an adiabatic 
plasma wave can be derived from the least action princi¬ 
ple, 5A = 0, where A is the action integral given by 


between 1C and the inertial frame, or remapping x 
merely introduces or removes a chirp determined by E.) 
Hence we will call /C the laboratory frame for simplicity. 
In case of a periodic system, /C can be the true laboratory 
system indeed (Appendix A). Otherwise, one can always 
choose 1C such that, locally, its velocity relative to the 
true laboratory frame be zero, even though the relative 
acceleration may be nonzero. Also keep in mind that 
the laboratory frame of reference is not necessarily the 
plasma rest frame. 

B. OC Hamiltonians 


A = 


J [(-Gem) - ' 


]ns{{ns)) 


dtdPx. 


( 1 ) 


Here £em is the Lagrangian density of electromagnetic 
field in vacuum, are the oscillation-center (OC) 
Hamiltonians of single particles of type s, the summation 
is taken over all particle types, hs are the corresponding 
average densities, (...) denotes averaging over rapid os¬ 
cillations in time and space, and ((...)) denotes averaging 
over the local distributions of the particle OC canonical 
momenta. Keep in mind that these distributions must 
be treated as prescribed when the least action principle 
is applied to derive the NDR. 

We will limit our consideration to electrostatic oscilla¬ 
tions in one-dimensional nonrelativistic electron plasma. 
In this case, £em = E'^/Stt, where E = —8x4) is the 
electric field, and </> is the potential. We will allow 4> to 
consist of a rapidly oscillating potential 4> and, possibly, a 
potential 0 that is “low-frequency” (LF) both in time and 
space. Depending on a specific problem (Appendix A), 
the LF field, E = —dx4>, may be zero, determined ex¬ 
ternally, or result from wave dynamics, in which case it 
can also be calculated, if needed. However, finding E is 
a problem that is separate from finding the NDR. This 
is because, strictly speaking, a local NDR is well defined 
only in the (possibly noninertial) frame of reference K. 
where the A-driven acceleration vanishes, at least locally 
[41]. Denoting the new coordinate as x, one can write 


A = 


£ J'[A] dtd^x. 


( 2 ) 


Here, the wave Lagrangian density, £, is given by 


(dxE)^ {{8x4)?) 
Stt Stt 




( 3 ) 


and J\E] is the Jacobian of the coordinate transforma¬ 
tion. Since the mapping a; i—a: is independent of 4), the 
fact that J is generally a functional of E has no effect on 
the NDR [42]. We can also omit the term {8x4 )Y/{^t^) 
Eq. (3) for the same reason. 

It is thereby seen that whether the LF field is zero 
or not is irrelevant to further calculations. (Switching 


The OC Hamiltonians of passing particles (denoted 
with index p) and trapped particles (denoted with index 
t), which enter Eq. (3), can be expressed as follows: 

Up = Pu + S — mu^/2, Tit = £ — mu^/2. (4) 

Here u = uj/k is the wave phase velocity (we use the sym¬ 
bol = for definitions), w = —8t4. is the local frequency, 
k = 8x4, is the local wave number, and 4 is the wave 
phase. Also, P is the OC canonical momentum of a pass¬ 
ing particle in the laboratory frame 1C, 

£ = mw^/2 + e4> (5) 


is the particle total energy in the frame 1C' that moves 
with respect to 1C at velocity u, w is the particle velocity 
in 1C', and m and e are the electron mass and charge. 
Notice also, for future references, that, at small enough 4), 
Up becomes the well known ponderomotive Hamiltonian, 


Up 


p2 

2m 


4>. 


( 6 ) 


Here the “ponderomotive potential” $ is given by 


$ = 


e^k'^4)\ 

4m{uj — kP/mY ’ 


( 7 ) 


and 4)1 is the amplitude of 4)- The OC canonical momen¬ 
tum P can then be expressed as 


P « mV — 9y4>, 


( 8 ) 


where V is the average velocity, or the OC velocity. 
Using Eqs. (4), one can cast Eq. (3) as 


£ = 


Stt 


n{{£)) - np{{P))u 


nmvf 

2 


(9) 


where fi = Up + fit is the total average density. Let us 
also introduce J as the particle action variable in K.', 
i.e., the (appropriately normalized) phase space area that 
is swept by the particle trajectory on a single period, 
J ex m § wdx (Fig. 1). The separatrix action, J*, can be 
estimated as J* ~ Jt/afi, where we introduced 


J = 


k^ ’ 


Op = 


ek'^f’i 

mujp 


( 10 ) 
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FIG. 1: Schematic of particle trajectories in phase space, illus¬ 
trating the definition of 2ttJ (shaded area): (a) for a passing 
particle, (b) for a trapped particle. For J to be continnous 
at the separatrix (dashed), with the passing-particle action 
defined as 2nJ = m§wdx, for trapped particles one must 
use the definition 27rJ = {jnl2) f wdx. The fignre is an ad- 
jnstment of Fig. 1 from Ref. [36]. 


property.) Apart from the examples considered in this 
paper, see also Refs. [28, 37, 45] where such decoupling 
is utilized to arrive at results of practical interest. 

Below we will present several models for F{J) and 
demonstrate that they readily yield the many scalings 
for 6uj that were reported previously for various special 
cases. We postpone considering beam distributions until 
Sec. V, whereas the bulk plasma will be modeled as fol¬ 
lows. Let us assume that a wave is excited by a driver 
with homogeneous amplitude, so the OC canonical mo¬ 
menta P are conserved and equal the initial momenta, 

P = mVo. (14) 


From Eq. (11), one gets 


Trapped particles, for which J is the OC canonical mo¬ 
mentum, have 0 < J < J*. Passing particles, for which 

P = mu + kJ sgn{w), (11) 

have J > J^. One can then rewrite Eq. (9) as follows: 

£ = ^ - n{{E)) + AH. (12) 

OTT 


J=\l-Vo/u\J, 
so the action distribution is given by 


F{J) = - 
m 


kJ\ 


fo[u + — j + fo[u - — 


kJ 


(15) 


(16) 


where /o is the initial velocity distribution. This is also 
known as the adiabatic-excitation model [10, 29]. 


Here AH is independent of the wave amplitude, so it has 

no effect on the NDR, as will become clear shortly [42]. Nonlinear Doppler shift 


C. Action distribution 

The ensemble averaging can be expressed through the 
action distribution F{J) that includes both trapped and 
passing particles, F{J) = F^{J) -|- Fp{J), where Ft_p(J ^ 
J*) = 0. We will assume the normalization such that 
/o°° dJ = fii^pln- then, F{J) dJ = 1, and 

pOO 

((...))=/ (...)F(J)dJ. (13) 

Note that the action distribution is not different from, 
say, the velocity distribution in the sense that, if needed, 
F{J) can be calculated self-consistently [43]. Yet F{J) is 
also advantageous in the sense that it is easier to calcu¬ 
late or model in many cases of interest, especially those 
that are of interest to approach analytically. This is be¬ 
cause J is the natural variable for phase-mixed distri¬ 
butions. For example, J is an adiabatic invariant for 
trapped particles, and passing particles often conserve P 
or T-Lp, which are expressible through J. Models of the 
distribution functions in such variables are therefore par¬ 
ticularly robust. Hence, the problem of modeling F{J) 
often can be decoupled from the problem of finding the 
NDR. In other words, it can be meaningful to search for 
uj for a given F{J). (The situation is similar to calculat¬ 
ing linear dispersion for a given “unperturbed” velocity 
distribution [44, Chap. 8], but nonlinear calculations in¬ 
volving velocity distributions usually do not enjoy this 


Before we proceed to a general calculation, notice the 
following. The conservation of P implies that, as a re¬ 
sult of the wave excitation, the plasma generally changes 
its average velocity, {{V)), by some {{AV)). This leads 
to a nonlinear Doppler shift, k{{AV)), by which the fre¬ 
quency UJ in the laboratory frame /C differs from that in 
the plasma rest frame [5, 46]. The shift can be calculated 
using that Vt = u and Vp = OpPp = u + {dj£)sgn{w), 
where £ is considered as a function of J and of the wave 
parameters [33]. However, the concept of the plasma rest 
frame is meaningful mainly when the resonant population 
is negligible, in which case the following simple deriva¬ 
tion is possible. (Also see Appendix B for an alternative 
derivation under the same condition.) 

By combining Eq. (8) with Eq. (14), one gets for an 
initially-resting plasma that 

((AE)) = {{V - Eo)) « ^ ^ (((1 - Eo/u)-3)). 

(17) 


One can also reexpress the right hand side as follows: 


m 


e'^k'^(t>\ 

Am^ 



kVo)-^]fo{Vo)dVo 


e^k4>l d /o(Eo) 

4m2 dujJ_^{Vo-uj/k)^ 
a£‘ uj^ de(uj, k) 

4 kujp duj 


dVo 


(18) 
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where 


E. Dielectric function 


e{uj,k) = l--^ 


/o(Vb) 


Loo {v^-^/kf 


dVn 


fc2 


mo) 

Vq — uj/k 


dVo. (19) 


This leads to 


((AF)) 


de{uj, k) 
4 ku!p du) 


( 20 ) 


where, within the adopted accuracy, the right hand side 
must be evaluated at the linear frequency, wq. Hence, 
one eventually arrives at the following expression for the 
nonlinear Doppler shift: 


km)) 


^ ^ dejup, k) 
4 a;2 (9wo 


( 21 ) 


where we introduced oq = ek'^(j)i/{mujQ). Also notably, 
the momentum density associated with the nonlinear 
Doppler shift is 


mn{{AV)) 


kE'^ de{ujQ, k) 
IOtt i9wo 


Wo ’ 


( 22 ) 


where E = k((>i is the wave amplitude, and € is the wave 
energy density [44, Sec. 4-4]. One may recognize k€/uJo 
as the canonical [47] momentum density of a linear wave 
[44, Sec. 16-3]. It is seen then that the nonlinear Doppler 
shift is precisely the effect that allows a linear electro¬ 
static wave carry momentum. (The additional accelera¬ 
tion caused by a LF field, if any, changes the momentum 
of the OC motion rather than the wave momentum.) 

Some calculations reported in literature ignore the 
nonlinear Doppler shift and present the NDR as a relation 
between Slo and gq in the plasma rest frame [48]. Keep 
in mind, however, that this frame is by default unknown. 
Hence, calculating Slo in the plasma rest frame gener¬ 
ally says nothing about the frequency shift that could be 
directly measured in laboratory. (Defining the wave mo¬ 
mentum is also problematic in this case.) Although it is 
in principle possible to excite a wave without accelerat¬ 
ing plasma, such scenarios are not typical for traveling 
waves that we consider [5]. Of course, there are geome¬ 
tries where the average flow is bound to vanish eventually, 
e.g., to prevent the growth of the LF field (Appendix A). 
But this effect is not universal and, in any case, is re¬ 
stricted to stationary waves. In contrast, our Lagrangian 
formulation allows calculating a more fundamental, local 
NDR, since ^ is treated as an independent function. As 
a result, our theory automatically accounts for plasma 
acceleration associated with the wave excitation in the 
frame K. for a given particle distribution. (As we ex¬ 
plained in Sec. HA, the possibly nonzero simultaneous 
acceleration caused by 4> is not a part of the NDR and 
can be found separately, if needed.) Hence, Sto that we 
calculate below already includes the nonlinear Doppler 
shift, and Eq. (21) will be used only to establish connec¬ 
tions between our results and other theories. 


It is to be noted that integration by parts in Eq. (19) 
is justified only by the assumption of having no resonant 
particles, in which case the integrand is analytic. More 
generally, we will define e as the following real function, 

eiLo,k) = 1-^P r J^dw, (23) 
k^ J_^v-uj/k 

where P denotes the Cauchy principal value. [In the 
absence of resonant particles, Eq. (23) is equivalent to 
Eq. (19).] For the adiabatic waves of interest, which are 
by definition phase-mixed and exhibit no Landau damp¬ 
ing, such e can be recognized as the linear dielectric func¬ 
tion [44, Chap. 8]. 


III. WAVE MODEL 
A. General dispersion relation 

Since large amplitudes typically result in rapid deteri¬ 
oration of waves through various nonlinear instabilities, 
the very concept of the NDR is of interest mainly for 
waves with small enough amplitudes, i.e., when the wave 
shape is close to sinusoidal. Harmonics of order £ > 1 
(which are phase-locked to the fundamental harmonic, 
£ = 1) can then be treated as perturbations, and the 
spectrum decreases rapidly with £ [49]. Typical calcula¬ 
tions of the trapped-particle nonlinearities, such as those 
reviewed in Ref. [36], are restricted to the sinusoidal wave 
model, which neglects all harmonics with £ > \. Below, 
we extend those calculations by introducing a “bisinu- 
soidal” wave model, which retains the first and second 
harmonic, while neglecting those with £> 2. 

Within the bisinusoidal wave model, one can search for 
the wave electrostatic potential in the form 

(^ = (^icos(c)-f 02COs(2^-f x), (24) 

where (pi, p 2 , and x are slow functions of {t,x). (We 
will assume that (p 2 is small enough, so there is only 
one minimum of the potential energy per wavelength; see 
below.) It is then convenient to adopt L X: and 

ai = epik^ / {mto‘^), 02 = ep2k‘^ / ifnio'^) (25) 

as four independent variables. Minimizing A with respect 
to ^ leads to a dynamic equation for the wave amplitude 
[33], which is not of interest in the context of this paper. 
The remaining Euler-Lagrange equations are as follows: 

da,^ = 0, da,^ = 0, d^£2 = 0, (26) 

where £ = £(ai, 02 , X) a;, k). By substituting 
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one hence obtains the following three equations, 


0 Uf» 



(28) 

(9a 1 

Sirn 

\ e / ’ 

n _^«^)) 

da2 

27rn 

\ e / ’ 

(29) 

n _^«^)) 

dx 



(30) 


where £ = J, oi, 02 , X; w, fc). After eliminating 02 and 
X and introducing a = ai, one arrives at a single equation 
for w(fc, a), which constitutes the NDR. 


B. Eliminating x 


Our first step is to eliminate x- To do this, consider 
Eq. (30) in the following form: 

((5xf(J,ai,a2,x,w,A:))) = 0. (31) 


In order to express £ as a, function of J, we will need an 
explicit formula for J. A universal formula that applies 
to both trapped and passing particles is [33] 


/•27r/fe ^_ 1 

J = Re / ^2m[£ — e<j)(x)\ —. 

Ja 27 r 


Using the notation 

^02 ^1 

2ai' 2 V rnu^a 


+ 1 


(32) 


(33) 


it is convenient to cast Eq. (32) as 


J = Jy/a J (^0/2) — z cos{26 + d9. 

Hence we can rewrite Eq. (31) as follows [50]: 

{{d^J/drJ)) = 0, (34) 

where J = J{r, a, z, x, w, k). Notice now that 


idxJ)\x=o 



dO 


X zsin(20) r — sin^(0/2) — 2 ;cos(20)] 


which is zero, because the integrand is an odd function of 
9. Therefore, Eq. (34) has an obvious solution, x = 0 (or 
X = TT, but the difference between these solutions can be 
absorbed in the sign of z). We will assume this solution 
without searching for others, because it corresponds to 
what seems to be the only (relatively) stable equilibrium 
seen in simulations, even for large-amplitude waves [51]. 
For an alternative argument, see [52]. 

The remaining equations that constitute the NDR 
[Eqs. (28) and (29)] can then be represented as follows, 


Qi 


aiuP' 


= 0, 


G 2 


2a20j‘^ 


UJt. 


= 0 , 


(35) 


or, equivalently, 


1 — 2Gilcip — 0 , 8zGi — G 2 — 0 . ( 36 ) 


Here Wp = (47rne^/m)^/^, Gi ,2 = ((G 1 . 2 )), and 


G 


1.2 — 


d 

dai^2 


£{J,ai,a 2 ,uj,k) 

mu^ 


(37) 


Notably, 9a2G'i(T,oi, 02 ,w, fc) = 9aiG2(T, ai, 02 ,w, fc) by 
definition, and the same applies to Qi^ 2 - Also notably, 
Gi _2 have a clear physical meaning, which is as follows. 
The function 2Gi/a can be understood as the relative 
contribution of a particle with given action (or energy) 
to the wave squared frequency, in units Wp. The 
function G 2/2 can be understood as the particle relative 
contribution to the second harmonic amplitude, (j) 2 , in 
units mojp/{ek‘^) [because a 2 W^/Wp = e 4 > 2 k‘^/{mojp)]. 


C. Normalized action 


It is convenient to represent the particle action as J = 
Jx/a j{r, z), where the “normalized action”, j{r,z), is a 
dimensionless function given by 


j{r,z) = -Re[ x/n{r,z,9) d9, 

^ Jo 

TZ{r,z,9) = r — sin^(0/2) — zcos{29). 

Then, G 1.2 can be expressed through j alone, 

Gi = 2r - 1 -I- 2ai 9air( J, Oi, 02 , w, k) 
daiJ{r, ai,a2,u],k) 


= 2r — 1 — 2ai 


= 2r - 1 - 


drJ{r, ai,a2,uj,k) 
j - 2zjz 
Jr 


G 2 = 2aida2r{J,ai,a2,u},k) 

002 J(r,ai,a 2 ,w,fc) 


— —2ai 

= 

• 5 

Jr 


drJ{r^ai,a2^oo,k) 


(38) 

(39) 


(40) 


(41) 


where j = j{r, z), and the indexes r and 2 ; denote partial 
derivatives. Explicit formulas for j{r,z) and Gi_ 2 (''’)^) 
are given in Appendix C. Characteristic plots of Gi _2 as 
functions of r are presented in Fig. 2. 


D. Restrictions on r and z: 
reflection points and the trapping condition 


Let us also introduce an alternative, all-real represen¬ 
tation of j. 


jir,z) = -[ x/n{r,z,9) d9. 

TT Jo 


(42) 
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FIG. 2: (Color online) Gi(r,z) (solid) and G 2 {r,z) (dashed) vs r at fixed a: (a) a = —1/8, (b) z = 0, (c) z = +1/8. Notice the 
singularity at r = 1 + 2 , which corresponds to the boundary between passing and trapped trajectories (Sec. HID). 


Here, for passing particles 0q equals tt and for trapped 
particles Oq equals the first positive solution of 

7^(r,2,0o) =0. (43) 

Equation (43) can be expressed as the following quadratic 
equation for y = cos Oq , 

Azy'^ — j/ — 2r + 1 — 22 = 0. (44) 

That gives y = (1 ± D)/{Sz), where 

D = v'l - 162 + 32r2 + 322^. (45) 


of 2 . Assuming 2 <C 1 (which is verified a posteriori), 
one can hence replace j(r, 2 ) with its Taylor expansion, 

00 

j{r,z) = '£z-Ur). (50) 

71 — 0 

The coefficients /„ are calculated (Appendix C) using 

2 i—Ij" d" /- 

^ ^ /o (20) - sin^ (0/2) dO. 

By also expanding these functions in l/r, one gets 


The solution for y must be combined with the conditions 

-l^y^l, -1/8^ 2 <1/8, (46) 

which flow, respectively, from the definition of y and from 
the assumption of having a single minimum of the po¬ 
tential energy per wavelength (Fig. 3), which we adopted 
earlier. That leaves only one of the roots, 


j(r, 2 ) = 


2 


1 3 + 42^ 5 + 32 + 12z^ 

^ 32^^2 I28r3 

5 (35 + 482 + 1442^) 

8192r4 ■■■ ^ 


(51) 


where we omitted terms of higher orders (in both 1 jr and 
2 ) as insignificant for our purposes. The inverse series is 


00 = arccos > (47) 

and leads to the following trapping condition, 

2 < r < 1 + 2 , (48) 

Notably, this condition also guarantees that D and 0o are 
real for trapped particles. 

E. Asymptotics 

Below, we will also need asymptotics of j and G' 1 , 2 . 
For deeply trapped particles {r ^ z), one can show that 

j(r-)> 2 ) = 0, G'i, 2 (r+► 2 ) =+1, (49) 

whereas, for passing particles far from the resonance 
(r 1), the asymptotics are derived as follows [53]. 
Away from the separatrix, j (r, 2 ) is an analytic function 



_-t_A_t_A_ij 

0 7r/2 TT 37r/2 2 tt 


FIG. 3: (Color online) Potential energy, e(j>, vs the wave phase, 
for 2 = 0 (dotted blue), | 2 | = 1/8 (solid red), and | 2 | = 1/3 
(dashed black): (a) 2 ^ 0, (b) 2 < 0. A second (per period) 
minimum of ecp appears at \z\ > 1/8. 
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then found to be 

'■■2 1 1 + 4^2 3z 5(1 + 32z2) 


’■<^■•*■> = ¥ + 5 ' sp ■ 8j 


;4 


128j 


;6 


This leads to the following equation for S, 
S {J / JY a\ + a\ 3afa2 


mu^ 


4(J/J)2 8(J/J)4 

5(a| + 8afa|) 




so Eqs. (37) lead to 


1 3z 5(1 + 16z2) 

~ 2f ^ 2j4 + 16 j6 

z 3 5z 


+ ... 

(52) 


(53) 

(54) 

(55) 


Note also that, to the extent that ¥ can be neglected, 
Eq. (52) can be cast in the following transparent form, 


iP' 


+ 


e^k‘^4>i 


2m 4m(A:P'/m)2 


(56) 


where we kept only the terms that will be relevant for our 
discussion (the applicability conditions will be discussed 
below; also see Appendix D) and introduced 


pOO 

Vi= [4'i(j)^)j - l/2]c(?, 

^0 

(59) 

pOO 

V 2 = ['i> 2 {j,z)j - z]dj, 

Jo 

(60) 


where z) = — fg G'i, 2 (¥(J), z) dj. Hence, the NDR 

[Eqs. (36)] can be cast as follows: 


0 = e(w, fc) + w 


2 d‘^e{uj, k) za 


■ (jj 


doj^ 2 
d^e{u},k) 
duj^ 


192 ~ 4771^0M^/o(u) (61) 


0 = 3z + w2^-^^^-2772V^uVo'(w)^- (62) 

The linear frequency is then found as a solution of 

e(a;o, k) = 0. (63) 


Here P' = —kJ, which serves as the canonical OC mo¬ 
mentum in 1C'. [This is seen if one writes Eq. (11) in 1C', 
where the frequency is zero.] Hence, the second term on 
the right-hand side is recognized as the ponderomotive 
potential (7) produced by a zero-frequency wave in /C', 
and the above formula for £ is recognized as an approxi¬ 
mate formula for the canonical OC energy, Eq. (6). 


Also one finds that z = o(l), namely, 
2 d^e{uJo, k) ao 2 


a^S 48 ' 


(where uq = LOo/k), and the nonlinear nonlinear fre¬ 
quency shift, Jw = w — Wo, is given by 


IV. SMOOTH DISTRIBUTIONS 
A. Basic equations 

First, let us suppose that F is smooth compared to 
G, i.e., has a characteristic scale much larger than J* 
(Fig. 4a), at least near the resonance. This allows ap¬ 
proximating 0i^2 as shown in Appendix D. Specifically, 
using Eq. (D23) with coefficients Cq inferred by compar¬ 
ing Eq. (D2) with Eqs. (54) and (55), one obtains 




'de{uJo,k) 

-1 r 

dujo 



- Wn 


d'^e(uJo,k) zao 




— OJ, 


d'^€(uJo,k) al 
° duj^ 192 


4?7iv^/o (uo) 




k^ 


The coefficients 771^2 can be evaluated at z = 0: 

poo 

/ ['ki (j,0)j-1/2] -0.27, 

Jo 


m 


(65) 


( 66 ) 


Gi ~ y<i [1 -e(w,fc)] -( 


d‘^e{uj,k) za 

aa ;2 Y 


d'^€{u}, k) d 


doj^ 


^-h4?7iv^MVo(w) y|, (57) 


G 2 ~ ap\ [1 — e(w, k)]z — oj 


2 9^e(w, k) a 


doj^ 16 


+ 2Tj2Vau^fo{u) ^ p 


(58) 


which is taken from Ref. [36], and 

poc 1 poo 

772 « / 4'2(j,0)jdj = - / G2ir{j,0),0)f dj 

Jo ^ Jo 

G2ir,0)j‘^{r,0) jrir,0) dr « 0.11, (67) 

which is calculated numerically. Below, we demonstrate 
how these equations are applied to specific distributions 
and show how our results correspond to those that were 
reported in the literature previously. 





FIG. 4: Schematics of the specific phase space distributions that are discussed in the present paper: (a) arbitrary smooth 
distribution; (b) cold limit; (c) finite-temperature waterbag (flat) distribution; (d) deeply trapped particles; (e) flat distribution 
of trapped particles. In the cases (d) and (e), details of the passing particle distribution are inessential as long as the beam 
nonlinearity (Sec. V) dominates. The dashed line is the boundary (separatrix) separating passing trajectories (outside) from 
trapped trajectories (inside). Also, x is the coordinate, \ = 2'n jk, and u = w -|- u is the velocity in the laboratory frame, K,. 


B. Fluid limit 


Then Eq. (63) gives wq = Wp, and Eqs. (68) and (69) give 


Let us start with the fluid limit, when kXjj is relatively 
small (here \d = VTjtOp is the Debye length, and vt is 
the thermal speed), while a is relatively large (although 
a <C 1 is still assumed), such that one can drop the half¬ 
integer powers of oq in the above formulas. In particular, 
Eq. (64) then yields 


2 3^e(wo, k) oo 
dujl 48' 


( 68 ) 


By substituting this into Eq. (65) and again ignoring the 
term Oq^^, we obtain 


5ijj 


^ ' de{uJo,ky ^ 

96 dujQ 


X 



d'^e{uJo,k) 


^ d^e{uJo,k) \ 
2 dojQ j 


(69) 


z = ao/8, Suj = U}pall2. (71) 

This result is in agreement with Ref. [5], but it may 
seem to be at variance with the widely known theorem 
[48] that a Langmuir wave in cold plasma exhibits no non¬ 
linear frequency shift whatsoever (in fact, even at large 
ao). The discrepancy is due to the fact that the men¬ 
tioned theorem applies to the plasma rest frame only. 
Our calculation, in contrast, is done in the laboratory 
frame, where the excitation of a plasma wave is accom¬ 
panied by plasma acceleration. That results in the non¬ 
linear nonlinear Doppler shift discussed in Sec. IID. One 
can easily check, by substituting Eq. (70) into Eq. (21), 
that Suj described by Eq. (71) entirely consists of this 
shift. Hence, recalculating Eq. (71) to the plasma rest 
frame predicts zero frequency shift, as anticipated. 


2. Waterbag distribution 


The effect 6uj oc is called a fluid nonlinearity. Note, 
however, that trapped particles also contribute to this 
nonlinearity in the general case, contrary to a common 
misconception. This is because, unless /o is vanishingly 
small near the resonance, removing trapped particles 
would leave a phase space hole, rendering /o abrupt and 
the above results inapplicable. (The corresponding mod¬ 
ification of the NDR can be quite strong and is discussed 
in Sec. V.) What makes the plasma act as a fluid then 
is not necessarily the absence of trapped particles but 
rather the smoothness of /o near the resonance. 

To our knowledge, Eq. (69) is the first generalization 
of fluid nonlinearities to plasmas with arbitrary e. Below 
we consider several special cases for illustration. 


1. Cold limit 


Now suppose the waterbag distribution discussed in 
Ref. [31]. Specifically, assume that the initial velocities 
of particles are distributed homogeneously within the in¬ 
terval (—v,v), where v is some constant (Fig. 4c). The 
corresponding dielectric function is calculated by directly 
taking the integral in Eq. (19) [54] and is given by 


e(uj,k) 


1 - 


^2 _ J^2y2 


(72) 


Then Eq. (63) yields ujq = ujp + k'^v'^, and Eqs. (68) and 
(69) lead to 


ao(3 -I- a) 


24(1-d)2’ 


Suj « woUq 


2 (6 -I- 9 a -I- w ) 


12(1 - d)3 


(73) 

(74) 


In the limit of cold plasma that is initially at rest 
(Fig. 4b), one has 


(70) 


where a = {kv/ujof. At d —>• 0, Eq. (74) reproduces the 
result that we presented in Sec. IV B 1 for cold plasma. 

One may notice that Eq. (74) is at variance with the 
result obtained in Ref. [31]. This is due to the fact that. 


t{uj,k) = 1 
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in Ref. [31], the NDR is calculated not in the laboratory 
frame but rather in the frame where the average velocity 
is zero; hence the nonlinear Doppler shift is not included. 
For the distribution in question, Eq. (21) gives 


in electron plasma with the Maxwellian distribution, 


k{{AV)) 


UJqO? 

2(1-d)2- 


Subtracting this from Eq. (74) leads to 5u} = woa^Q;(15 + 
d)/[12(l — d)^], which reproduces the result of Ref. [31]. 


3. Kappa distribution 


We also calculated the effect of the fluid nonlinearity 
on the dispersion of a nonlinear Langmuir wave [55] in 
electron plasma with the kappa distribution [56], 


foiVo) 


r(^ + i)/r(^-i/2) r i " 

vtK\/ 7t{2k — 3) (2k — 3)w|. 

(75) 


where k > 3/2 is a constant dimensionless parameter, 
v^/o(w) dv is the thermal speed squared, and F 
is the gamma function. [The distribution (75) approaches 
the Maxwellian distribution at k 1, but has more pro¬ 
nounced tails at smaller k.] The results are shown in 
Fig. 5. Notice that 2 : is not monotonic in k, and Suj can 
have either sign. This is in contrast to Ref. [31], where 
the fluid frequency shift was reported as strictly positive. 


In this case, 

where Z{Q = —2S{(), and S is the Dawson function, 
S'(C) = exp(-C^) exp(j/2) dy [44, Chap. 8]. 

Two sets of results are presented. Figure 6 shows the 
results obtained by application of the asymptotic formu¬ 
las (64) and (65). Figure 7 shows the result of the direct 
numerical solution of the more precise NDR, Eq. (36). 
It is seen that the asymptotic formulas provide a rea¬ 
sonably accurate approximation to the exact solution, 
even though the applicability conditions of the asymp¬ 
totic theory, (D24), are satisfied only marginally. In par¬ 
ticular, the following trends are seen. At small kXo, the 
fluid nonlinearity dominates, leading to Eqs. (71). (The 
fact that, in the figure, Slo remains finite at vanishing 
kXo oc vt is due to the choice of the units in which 
the amplitude is measured.) At larger kXo, the kinetic 
nonlinearity becomes dominating. That said, compari¬ 
son with Fig. 5 shows that the presence of the kinetic 
nonlinearity does not affect the picture qualitatively. 


V. BEAM DISTRIBUTIONS 


C. Kinetic limit 


Now let us consider the opposite limit, when the am¬ 
plitude is relatively small, whereas kXo is substantial. In 
this case. 


z = 


2 fo (uo) 




Su « 4?7iv^/o (uo) 


woWp de{uJo, k) ^ 
kd' _ dujQ 


(76) 

(77) 


In addition to fluid and kinetic nonlinearities discussed 
above, waves with trapped particles can also exhibit non- 
linearities of a third, beam type. Those emerge when, 
on top of a distribution Fq that is smooth or negligi¬ 
bly small near the resonance, an additional distribution 
Fb of a trapped beam is superimposed that has abrupt 
boundaries within or at the edge of the trapping island. 
Such distributions, with both positive Fb (clumps) and 
negative Fb (holes), are known in various contexts [22- 
28, 37, 45, 57]. Below, we consider some of them in detail. 


The effect Slo oc y/oo is called a kinetic nonlinearity 
(associated with a smooth distributions; for abrupt dis¬ 
tributions, see Sec. V). In particular, Eq. (76) shows 
that (()2 oc This is in qualitative agreement with 

the (not self-consistent) estimate in Ref. [12], but notice 
that our numerical coefficient is different. As regard¬ 
ing Eq. (77), it is in precise agreement (modulo typos) 
with Refs. [34, 36] and also with the results reported in 
Ref. [10] for the “adiabatic excitation”. 

D. Fluid and kinetic nonlinearities combined 

We also calculated the combined effect of the fluid and 
kinetic nonlinearities for a nonlinear Langmuir wave [55] 


A. Basic equations 

The presence of Fb results in the appearance of addi¬ 
tional terms, 

Gi,2 = J GiM J) dJ, (80) 

in Eqs. (57) and (58). For simplicity, we will suppose 
that Fb is large enough, such that its nonlinear effect on 
the NDR dominates over the nonlinear effect produced 
by Fq. Then, 

C/i « (ap/2 )[l-e(w, fc)]- 1 -^ 1 , (81) 

^2 ^ apz[l - e{LO,k)] + G 2 , (82) 
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FIG. 5: (Color online) The effect of the fluid nonlinearity [Eqs. (68) and (69)] on the dispersion of a nonlinear Langmuir wave 
in electron plasma with the kappa distribution, Eq. (75): (a) normalized amplitude of the second harmonic, 2 : = 02/(201), in 
units Up; (b) nonlinear frequency shift, Slo, in units cupOp. The horizontal axes are kXo and the distribution parameter k. At 
K S> 1, the kappa distribution is approximately Maxwellian. 



FIG. 6: (Color online) The combined effect of the fluid and kinetic nonlinearities on the dispersion of a nonlinear Langmuir wave 
in electron plasma with the Maxwellian distribution, Eq. (78): (a) normalized amplitude of the second harmonic, 2 = 02/(20i); 
(b) nonlinear frequency shift, 5uj, in units cup. The horizontal axes are kXo and eExj[rnijjfVT) = Uosc/ut, where Ei = k4>i 
is the amplitude of the first harmonic of the wave electric field, and Uosc = eExj (mujp) is the characteristic amplitude of the 
electron velocity oscillations. The results are obtained by application of the approximate formulas (64) and (65). 


SO Eqs. (36) yield the NDR in the form 


B. Deeply trapped particles 


2 ~ 


Sop ’ 


5oj 


2^1 ' dejuJo, k) 

ttp [ dujQ 


(83) 


Suppose that trapped particles reside at the very bot¬ 
tom of the wave troughs (Fig. 4d); i.e., 


Notice that the wave is close to linear in this case not 
when Op is small, like before, but rather when /3 = 
fit,/{nap) is small, where fit, is the beam density. (In 
other words, Op must be sufficiently large.) The NDRs 
for specific distributions are derived as follows. 


Fb{J) = {nb/n)6{J). (84) 

[Strictly speaking, the delta function here must be un¬ 
derstood as the limit of 6{J — Jc) at —>■ 0.] Hence, all 
trapped particles have r = 2 , so Eqs. (49), (80), and (84) 
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kXn 

FIG. 7: (Color online) Same as in Fig. 6 but using the more 
that are specified in the legend. 



NDR, Eq. (36), and only for selected values of eE/(mujpVT) 


yield Qi ^2 = T^h/n. Then, Eqs. (83) give 


2. Discussion 


P Sui 2/3 de(uJo,k) 

3 ’ u!p Wp [ dujo 


(85) 


The obtained expression for Sw agrees with the well 
known result [22, 23]. It is also easy to see that the 
value of z is in agreement with what one can infer from 
the Vlasov-Poisson system directly [58]. 


Since Qi given by Eq. ( 86 ) is independent of w, one 
can also derive a more precise analytic expression for uj 
from Eq. (81) when e is a simple enough function. For 
instance, for cold plasma [Eq. (70)], one gets 


(jj 

LOp 


1 + 


8Af 

37r 


- 1/2 


= 1 - 


4Af 

Sir 


8J\f^ 

37r^ 


+ 0(A/'3). (91) 


C. Homogeneous clumps and holes 

1. Basic equations 


Let us also consider the case when F^, is flat across the 
whole trapping area (Fig. 4e). In this case, Ft,{J < J*) = 
F, where F is some constant. Then, 


1 + 2 

Gi^ 2 {r, z) jr{r, z)dr 

= apM-tia{z) pz apM^i^2{0), ( 86 ) 

where we introduced the dimensionless quantities 





k Da 
l+z 


/ 3 , 


v 


pL-\-Z 

71 , 2 ( 2 )=/ Gi^ 2 ir,z)jrir,z)dr. 

J Z 

Hence, one obtains 


z 


^72(0)A/', 


Sw ~ 271 (O)A/' 


de{uJo,ky ^ 

duJo 


(87) 

( 88 ) 


(89) 

(90) 


We find numerically that 72 ( 0 ) ~ —0.085, and 71 (0) = 
— 4 /( 37 r) was reported already in Ref. [36]. 


Notably, phase space clumps {Af > 0) in this case cor¬ 
respond to Jo; < 0 , whereas phase space holes {Af < 0 ) 
correspond to Jo; > 0 . 

Let us compare results predicted by Eq. (91) for w, 
as well as those by Eq. (89) for z, with the exact kinetic 
solution, which happens to exist in this special case. [The 
word “exact” here refers to the description of trapped 
particles, whereas the bulk plasma is still modeled with 
a linear dielectric function (70).] Assuming the notation 
a = {tt/2){uJp/uj), the exact solution is given by [59] 


TT^ — 4a^ 


327r2 - 8a2 ’ 


A7-2 = 


2 tan a 
3a 


or, more explicitly for small Af, 


-4a4’ 


^ = -^ + 0{Af\ 


U) 


= 1 - 


457r 
43V 68A72 

277r2 


-0{Af^). 


(92) 

(93) 


(94) 

(95) 


As seen in Fig. 8 , the results of our asymptotic the¬ 
ory are virtually indistinguishable from the predictions 
of Eqs. (92) and (93) at |A/’| < 1 and even beyond. The 
asymptotic theory also captures the fact that there is no 
solution for z and uj below some threshold, Af < Afc] 
i.e., deep enough hole modes cannot stable [60]. Our 
Eq. (91) predicts Afc = — 37 r /8 ~ —1.18. This dif¬ 
fers by less than 4% from the exact-theory prediction. 
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FIG. 8: (Color online) The effect of the beam nonlinearity on 
the dispersion of a nonlinear Langmuir wave in cold electron 
plasma for the case when the trapping islands contain homo¬ 
geneous clumps (A/” > 0, as in Fig. 4e) or holes (A/” < 0): (a) 
normalized amplitude of the second harmonic, z = (j> 2 l{ 2 (j>i); 
(b) nonlinear frequency shift, Sco, in units cOp. The coordinate 
Af, which is given by Eq. (87), represents a dimensionless 
measure of the trapped-particle phase space density. The red 
solid curves correspond to the exact solutions, Eqs. (92) and 
(93). The blue dashed curves correspond to our asymptotic 
solutions, Eqs. (89) and (91). There is no solution for z and 
a; at Af < Afc- The exact theory predicts A/), « —1.22, and 
the asymptotic theory predicts Afc ~ —1.18. 


wave second harmonic. It is shown for the first time that 
the fluid nonlinearity, which is determined by the second 
harmonic, can be cast in terms of the plasma dielectric 
function [Eq. (69)] and, contrary to the common pre¬ 
sumption, can have either sign. Also, the kinetic nonlin¬ 
earity that we derive accounts for both smooth distribu¬ 
tions [Eq. (77)] and trapped-particle beams (Sec. V). Our 
formulation is benchmarked against the many previously 
known NDRs and is shown to reproduce them as special 
cases of a single unifying theory. For example, we calcu¬ 
late the NDRs for Langmuir waves in plasmas with cold, 
waterbag, kappa, and Maxwellian distributions. These 
specific results and our general method are applicable, 
e.g., to waves produced at intense laser-plasma interac¬ 
tions and to energetic-particle modes in tokamaks. Our 
asymptotic formulation, in fact, may be advantageous 
over exact kinetic solutions for such waves. While offer¬ 
ing a reasonable precision, it leads to simpler (and thus 
more elucidating) results, does not involve solving any 
differential equations, and can be extended straightfor¬ 
wardly to other nonlinear plasma waves. 

One of the authors (lYD) thanks D. Benisti for valu¬ 
able discussions. The work was supported by the 
NNSA SSAA Program through DOE Research Grant No. 
DE274-FG52-08NA28553 and by the DOE contract DE- 
AG02-09GH11466. 


Me = —\/3/2 ft! —1.22, which is obtained from Eq. (93) 
in the limit a —>■ 0, i.e., at ui/wp —>■ oo. 

One conclusion from this is that the long-range sweep¬ 
ing of energetic-particle modes, which is described in 
Ref. [27] by means of the aforementioned exact solution, 
can be described with fidelity just using a simple asymp¬ 
totic NDR, Eq. (91). [In fact, Eq. (91) was proposed 
already in Ref. [26].] Using the theory presented here, 
asymptotic NDRs can be derived also when exact ana¬ 
lytic solutions do not exist. For example, in addition to 
beam nonlinearities, one may need to account for fluid 
nonlinearities when ^ M, and for kinetic nonlineari¬ 
ties when Up ^ AA^ and kXn 1. That is done simply 
by retaining all the terms in the expressions for Gi, 2 , in¬ 
cluding not only C/i 2 but also the other terms that appear 
in Eqs. (57) and (58). To the leading order, the nonlin¬ 
ear effects stemming from fluid, kinetic, and beam effects 
then enter the NDR additively. 


VI. CONCLUSIONS 

In summary, we derived a transparent asymptotic ex¬ 
pression for the nonlinear frequency shift of intense Lang¬ 
muir waves in general collisionless plasma. Our result 
describes kinetic and fluid effects simultaneously. This 
is achieved by application of a variational approach sim¬ 
ilar to the approach used in Refs. [33-36]. In contrast 
to the previous calculations, however, the present one 
takes into account corrections to the NDR caused by the 


Appendix A: Low-frequency field 

Here, we discuss some issues regarding calculations of 
the LF field, E = —dx4>, in electrostatic plasma approx¬ 
imation (EPA). A common way to define such an ap¬ 
proximation is to omit, in Ampere’s law, the curl of the 
magnetic field transverse to the current, Bj^. In projec¬ 
tion on the current axis, a:, this gives 

dtE + 4TTj = 0, (Al) 

where, unlike in the rest of the paper, j denotes the 
current density. For instance, this form was used in 
Refs. [6, 15]. It was concluded then, in particular, that 
dtE = —4:Tr{j), in which case a stationary wave is possible 
only at zero average current. 

However, strictly speaking, the assumption of one¬ 
dimensional dynamics implies that has a vanishing 
effect on the particle motion but not a vanishing curl. 
Those are not equivalent requirements. For instance, 
consider a cylindrical plasma column with current along 
the axis of symmetry. On the axis, one can have nonzero 
component of V x B^ along the current even though Bj^ 
itself is locally zero. Off the axis, Bj^ is a continuous 
function of the radius R, so its effect is small at small 
enough R, while the curl of Bj^ is nonnegligible. (This 
is a standard assumption, e.g., in the plasma probe the¬ 
ory.) Moreover, the effect of Bj^ can be made arbitrarily 
small everywhere, if an additional strong homogeneous 
magnetic field is imposed parallel to the plasma column. 
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This shows that Eq. (Al) is not always reliable for de¬ 
termining E in the EPA. Gauss’s law can be used instead, 

= -dTT/J, (A2) 

where p is the charge density. (Note that, when calculat¬ 
ing p on envelope scales, one generally needs to account 
for ponderomotive forces.) But Eq. (A2) contains spatial 
derivatives, so it must be supplemented with boundary 
conditions. These boundary conditions are determined 
by external charges at plasma walls (if any), so they 
are independent of the average current. This means, in 
particular, that the nonlinear Doppler shift described in 
Sec. IID may or may not be associated with the build-up 
of a LE electric field, depending on a particular system. 
Eor instance, absent external inductive fields, the aver¬ 
age field in case of periodic boundary conditions (e.g., in 
toroidal plasma) will remain zero at all times [61]. 

Note also that Eq. (Al) generally violates the funda¬ 
mental energy-momentum conservation in EPA, as dis¬ 
cussed in Ref. [62]. In contrast, the field-theoretical 
approach developed in the main text, which implies 
Eq. (A2) yet not Eq. (Al), predicts that the wave exci¬ 
tation inputs precisely the amount of canonical momen¬ 
tum expected from the linear wave theory (Sec. IID; also 
see Ref. [33] for the energy conservation). It is possi¬ 
ble that, eventually, a LE field is established and redis¬ 
tributes plasma momentum or transfers some of that to 
the walls; however, as we explain in Sec. II A, these effects 
are unimportant to calculating a local NDR. Given an 
OG distribution for a specific time of interest, the NDR 
can be calculated irrespective of the factors that have 
shaped that given distribution. In this sense, whether E 
is zero or not is irrelevant in the context of this paper. 

Appendix B: Bulk plasma acceleration by 
homogeneous wave excitation 

When a homogeneous electrostatic wave is excited in 
plasma, the plasma generally changes its average veloc¬ 
ity by some {{AV)). One derivation of this effect was 
presented in Sec. IID using the OG formalism. Here, we 
show how a hydrodynamic model yields the same result. 

Gonsider the equation for the plasma momentum den¬ 
sity, V = mnv, namely, 

dtV + dx{mnv^) = neE — (Bl) 

where n is the particle density, v is the flow velocity, and 
n is the pressure. By averaging over x, which we denote 
with (.. .)x, one gets 

dt{'P)x = {neE)^ « {heE)x- (B2) 

(The tilde in this appendix denotes quantities of the first 
order in the field amplitude.) Here we used that the gra¬ 
dient of a periodic function averages to zero. We also as¬ 
sumed zero external LE field and periodic boundary con¬ 
ditions, so a LE field remains zero at all t (Appendix A). 


Notice that eE that enters Eq. (B2) can be found from 


eE = mdtv + dxti/fi- (B3) 

Assuming H = n(n), we have n dxil/n = const x dxfi^, so 

{neE)x = m{ndtv)x = mdt{nv)x - m{vdtn)x- (B4) 

From the continuity equation, we also have dth = —ndxV- 
This leads to 

{heE)x = mdt{hv)x + mn{vdxv)x 

= dt{mnv)x + {dximnv'^/2))x 
= dt{mnv)x- (B5) 


Thus, dt{V)x = dt{mnv)x, or, after integrating over t, 

{AV)x = {mhv)x = mu{n^)x/n = mn{v^)x/u, (B6) 

where we also substituted the continuity equation, n « 
nv/u. Then, {{AV)) = {A'P)x/{mn) is given by 

{{AV)) = u{{n/n)‘^)x = {v^)x/u. (B7) 

Since, in the assumed model, one has e(a;, A:) = I — 
LOpKio'^ — 3k^v^), where 3v^ = dnn{n)/m [44, Ghap. 3], 
Eq. (20) is then easily recovered. 

Appendix C: Formulas for j and Gi ,2 

The explicit formulas for the auxiliary functions intro¬ 
duced in the main text can be expressed in terms of the 
following special functions, 

/.7r/2 

K{C)= (l-Csin2 6»)-i/2^6l, 

Jo 

p77/2 

E(C)=/ {l-Csin^e)^/^de, 

Jo 

^7r/2 

n(d,C)= / (l-dsin2 6l)-i(l-Csin2 6»)-i/2rf6l, 

Jo 

which are the complete elliptic integrals of the first, sec¬ 
ond, and third kind, respectively [53]. 

1. Normalized action 

The explicit formulas for the function j are as follows. 
For trapped particles {z ^ r < 1 + z), one has 


- KK(Rt) + /IfE( r,) + /i,"n(«„ R,)], 

(Cl) 

* 80 -b 1 - D ’ 

(G2) 


(C3) 

* Sz + l-D ' 

(C4) 

2r + 6z — 1 + D 

Jjvf — • 

2D 

(C5) 

2D 

(G6) 
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For passing particles {r > 1 + z), one has 


- [AKK(i?p) + AEE(i?p)-bAnn(iVp,Rp)] , 

(C7) 

— z — 1) 8z — 1 + D 

P V2r-b6z-l-bH 8z-bl-bT>’ 

(C8) 

E 12r + 6z — 1 + D 

p - V 2 

(C9) 

2y2(r-z-l) 1 

P V2r -b 6z - 1 -b G 8z -b 1 -b ’ 

(CIO) 

f? 

^ 2r -b 6z - 1 -b H ’ 

(Cll) 

P 2(r - z) ■ 

(C12) 


Note also that, for passing particles, D is real for all r 
when z ^ 0 and for r ^ 1/2 + \z\ + 1/(32|2;|) when z < 0. 
Even when D is imaginary, however, the above formulas 
still yield j as a real function. In particular, the first 
three coefficients in Eq. (50) are real functions given by 

/i = -V^E(r-i), (C13) 

TT 

h = ^ ^ ^ [16r(-l + 2r)E(r“^) 

-2(3-16r + 16r2)K(r"i)], (C14) 


107r(r — l)y/r 

X [(-5 + 224r - 1248r2 + 2048r3 - 1024r‘^)E{r-^) 

+ 16(5 - 47r + 138r2 - IGOr^ + 64r^)K(r-i)]. 

It is to be noted that, in the limit 2 —>■ 0, these formulas 
are commonly known. For example, see Ref. [36]. 


2. Functions Gi 2 


The explicit formulas for the functions Gi ^2 are as fol¬ 
lows. For trapped particles (z < r < 1 -b z), one has 


Gi = A]' 


inn(iVt,Rt) 

K(i?t) ’ 


A'E(R0 + A^n(fVt,Rt) 

Gr2 — -C>t H- 


4 = - 


K(i?t) 

1 + D 




8z ’ 

1 - 8z -b D 

8z ’ 


-K 1 - 32z 2-b (1-b 8z)i:) 

* “ 32z2 

* 2z 

l-Sz + D 
‘ 32z2 ■ 


(C15) 

(C16) 

(C17) 

(CIS) 

(C19) 

(C20) 

(C21) 


For passing particles (r > 1 -b z), one has 


G 2 = B'2 + 


Ann(iVp,Rp) 

+ K(Rp) ’ 


K(Rp) 

l-D 


= - 
P 8z 


^P = 


1 -b 8z - H 
8z ’ 


-K 1 - 8z -b 16rz -b 16z2 - D 
Bp = - 


B^ = 


32z2 

1 — 2r — 6z — H 


=- 


4z 

1 -b 8z - 


32z2 


(C22) 

(C23) 

(C24) 

(C25) 

(C26) 

(C27) 

(C28) 


It is to be noted that, in the limit z —>■ 0, the function 
Gi was reported also in Ref. [36]. 


Appendix D: Approximations for integrals over 
smooth distributions 


Here, we show how to approximate integrals like 

pOO 

g= GF{J)dJ, (Dl) 

Jo 

where F is smooth compared to G. We assume that G = 
G(r(j)) = G(j) is either Gi or G 2 , so it has asymptotics 
[cf. Eqs. (54) and (55)] 


G(j) 


C2 


3c4 

f 



(D2) 


where Cq are independent of j. (The dependence on z 
is implied but will not be emphasized in this appendix.) 
Then, an approximate formula for Q is derived as follows. 
First of all, let us introduce an auxiliary function 


Y{J) = - [ GdJ. (D3) 

Jo 


It is convenient to express it also as F = where 


4> = -fG{j)dj (D4) 

Jo 

is a dimensionless function with the following properties: 


^f(j ^ 0) = 0, 4'(j 00 ) = 0. (D5) 


The former equality is obvious from the definition, and 
the latter equality is understood as follows. Notice that, 
with G being either of Gi {i = 1,2), the function 'I'(j) 
is proportional to Gi calculated on the waterbag (flat) 
distribution with thermal speed v oc j. At v —>■ 00 , the 
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boundaries of this distribution are infinitely far from the 
resonance, so the trajectories on these boundaries are un¬ 
affected by the wave. Thus, {{S)) is insensitive to the wave 
amplitude, and thus both Qi and G 2 are zero. [For z = 0, 
this is also shown in Ref. [36] via a straightforward cal¬ 
culation.] Thus, in the limit of large j, the function 'l'(j) 
must vanish, which proves the second part of Eq. (D5) 
and also leads to the following asymptotics at large j: 

'l’U) = - + ^ + ^ + ... (D 6 ) 

J r 3 

Hence one can rewrite G as follows, 


approximate the integrand, Q{J) = Y{J)F'{J), as 

Q{J)=Qi{J) + Q2iJ)-Q3iJ), (D 8 ) 

where the indexes denote the small-J, large-J, and 
intermediate-J asymptotics, correspondingly. At small 
J, Qs cancels Q 2 , so Q(J) « Qi{J)- At large J, Q 3 can¬ 
cels Qi, so Q{J) « Q 2 {J)- Thus, Eq. (D 8 ) approximates 
the true function Q accurately at all J. 

To construct the specific expressions for notice 

that the derivatives of F at J = 0, = F(‘?^(0), are 

yielded by Eq. (16) in the form [63] 


G = - 


pOO pOO 

/ Y\J)F(J) dJ= Y(J)F\J) dJ. (D7) 
Jo Jo 




k\ 


9+1 


0 , q is odd, 
1 , <7 is even. 


(D9) 


Using the fact that F' is smooth compared to Y, we can Then, much like in Ref. [36], one can take 
_I 


gi (J) = r (J) (^Fg^^) J + Fg^") ^ +...), (DIO) 

Q2{J) = J^a(^j + ^J\+^ + ... ^ F'( J), (Dll) 

Qs{J) = ^ Pa + ^ PP + ... ^ (Fg(^^ J + F^^^^ + Fg^®) + ... ^ (D12) 

After substituting this into Eq. (D8) and rearranging the terms, we can express Q as 


Q = J^a 


F'{j) (7 + |r +<■+)l ++) - (i| 1 ++^ jv) - fyT (^ jv) 


7'(2) / 


. g ^(Po) [^{j)j - C2] + Fg^'‘^(Pa)^ ^ [^'(j)j^ - C2f - C4] + P5°^(Pa)^^ [^'(j)j® - C2P - C4P - cg] + ... 




Hence, we obtain 


G = [(TP)c 2X^'^ + (J'a)'c4xW + + ... ] 

+ [(Pa)3/2Fg(")r7(2) + +...], (D13) 

where the coefficients are defined as the following (converging) integrals: 


J-^F'iJ)dJ, (D14) 

^0 

poo 

p4)=/ J-3[F'(J)-Fg("P]dJ, (D15) 

^0 

pOO 

^(6) ^ _ p(.2)j _ Fg(")p/3,] 

^0 

pOO 

77 ( 2 )=/ [^'(j)j - C 2 ]dj, (D17) 

^0 
1 7°^ 

= J\ j i'^i3)f - C 2 f - C 4 ] dj, (D18) 

1 

17 ^®^ = ^ P - C 2 f - C 4 f - cg] dj. (D19) 
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It is convenient to express in terms of the dielectric function given by Eq. (23). This is done as follows: 

= — 


/■“ fl T 

J [fo{u +kJ/m) - fo{u-kJ/m)] — 

k'^ rhn 


m 

u 2 


^2 




dw 


V — u 




[1 - e(w,fc)]. 


^ 


Using the formulas derived in Appendix E [namely, Eqs. (E7) and (E12)], we also obtain 

POO 7 T 

■ J [foiu + kJ/m) - fgiu - kJ/m) - 2fl^‘^\u) (kJ/m)] — 

r»4 poo 7 

= iiSJ J + w)-foiu-w)- 2 f^'^\u)w] 

k'^ p \f'( , 'i r(2)/ 1 

= + w)-h («)H ^ 


- Ui A_p 

m"* 2 dv? 


SMdv 


_ -aoV-U 

1 fc® d‘^e{uj,k) 

2 m'^LUp duj"^ ’ 


(D20) 


(D21) 


i-2 

x''> = 4 

W JO 

= «■" 

TO° e-s-O 


J fo{u + kJ/m) - fo{u - kJ/m) - 2 f^^\u) kJ/m - 2/^(m) 


dJ 

J5 


/o(u + w) - /o(m - w) - 2f^^\u)w - 2f^^\u) ^ 


3! 






=ffp r 

7-00 

_ 1 ^4 

4! 

.12 


/o(u + w) - /y (w)w - /y (w) y- 


„31 


3! 


dw 


foiv) 


dv 


_ k^^ d^e{uj,k) 

4! m®Wp 9a;4 

Then one can rewrite Eq. (D13) as 

■2^ d^e{w,k) 


Q = 


[1 - e{w, k)]c2 - 


2! aa;2 ' 4! dw'^ 

+ 2a|r?(^)ai/2M3^ ,;(6)a5/2^7^ ^ _ 

I- 


i9'^e(a;, k) ^ 


cea 


(D22) 


(D23) 


Einally, notice the following. In Eq. (D23), the term 
in the first curly brackets, which we call the fluid term, 
contains an expansion in integer powers of a. If Cq with 
neighboring q are about the same, the ratio of the neigh¬ 
boring terms is of the order of a. The term in the second 
curly brackets, which we call the kinetic term, contains an 
expansion in half-integer powers of a. If with neigh¬ 


boring q are about the same, the ratio of the neighboring 
terms is about au^/v^ (assuming, e.g., the Maxwellian 
/o), where vt is the thermal velocity. Hence, the appli¬ 
cability conditions of Eq. (D23) are 

a <C I, au'^/v^ <C 1. (D24) 

That said, the applicability conditions may not actu- 
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ally be as strict as this formal estimate suggests. That is 
because, when the fluid term is much bigger than the ki¬ 
netic term, the latter does not need to be calculated with 
high accuracy, and vice versa. It is also to be noted that, 
when deriving the conditions (D24), we entirely neglected 
order-one numerical coefficients. In practice, these coef¬ 
ficients seem to make higher-order terms less important. 
We conclude this from the results reported in Sec. IV, 
where our asymptotic formulas demonstrate reasonable 
agreement with more precise calculation based on the 
original NDR, Eq. (36), even though the applicability 
conditions of the asymptotic theory are satisfied for the 
chosen parameters only marginally. 


where h is some function, and where the following nota¬ 
tion is introduced, 

p pU — £ pOO 

{...)dv= / {...)dv+ {...)dv. (E2) 

J £ J —OO J U-\-£ 

Note one special type of such integrals, 

(2^-1)/? - = (E3) 


Appendix E: Derivatives of principal value integrals , . , . • , r i ^ i i , 

which IS a convenient formula to be used below. 


In this appendix, we report some useful identities in¬ 
volving derivatives of principal-value integrals. It is likely 
that these formulas can be found in the existing litera¬ 
ture, but it is instructive to rederive them here. Specifi¬ 
cally, we will consider integrals of the type 


P 



V — U e-i-0 V — U 


(El) 


1. Second derivative 


Consider an expression of the form 


1 

= — hm 

2 dy £-10 


h{v) 


dv 




V — u 


f U+£ 


h{v) 


V — u 


dv 


= i A lim 

2 OU e-iO 


= lim 
£—^0 


= lim 
£—^0 


h{u — e) 

-£ ^ J-oc {v-uf 

h'(u — s) h{u — e) f h{v) 


dv 


hiv) , h{u + e) 
dv - 


-2e 

h{v) 


2^2 

-du -f AP2 


(v — u)' 


■ dv — 


u+e (v - uY e 

h{u + e) h'{u-\-e) 


2e2 


2e 


(u — uY 

Here we introduced 

h'{u — e) h{u — e) h{u-\-e) h'{u + e) 2 h'{u) 


AP, = - - 

2e 2£2 2e2 

where Eq. (E3) was substituted. Hence, 

h{v) 


P 2 \h] = lim 
£^0 


dv — 


2 e 


h'{u) 


+ 0(1) = - f , ^ ^ dv + 0(1), 

Je Y - up 


dv 


{v — uY Je Y ~ '^Y 

This leads to the main result of this subsection, which is that 


= P 


h{v) h'{u) 


{v — uY {v — uY 


dv. 


i^p r 

2 ! dy V — u 


dv = P 


h{v) — h'{u){v — u) 
{v — uY 


v. 


(E4) 


(E5) 


(E6) 


(E7) 
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2. Fourth derivative 


The fourth derivative is treated similarly: 

1 ^ f°° h{v) 


1 d'^ f 

2 92 


dv 


4 ! du'^ 


P2[h] 


= — — P I 

12 (9u2 J 
= J_ A lim 

12 du e-^0 


= 1^1“ I# 

12 £->o 1 du 


h{v) 

h'{u) 

_{v — uY 

{v — uY 

Sh{v) 

2h'{u) 

_{v — uY 

{v — uY 

Sh{v) 

h"iu) ' 

_{v — uY 

{v — uY _ 


dv. 


dv + AP4 


dv + Tr- AP 4 

ou 


(E8) 


[the term 2h'{u)(v — u) ^ was dropped because the principal value integral over it is zero], where 


AP4 = 


h{v) 

h'{u) 


h{v) 

h'{u) 

_(u — uY 

{v — uY. 

v—u—s 

_(v — uY 

{v — uY. 


= —^ [h{u — e) + h{u + s)]. (E 9 ) 


v—u-\-s 


Then, 


P4 [h] = lim I [ 
12e^o\duJ^ 

where we introduced 

Ai ?4 = 

It is easy to see that 


3 h{v) h"{u) 


{v — uY 

3 /i(u) h"(u) 

{v — ■u )4 (v — u)2 

1 


dv — AP4 > = lim 
ou I e ->0 


hjv) h"'{u) 

(v — u)^ 12(u — m)2 


dv 


AR4 


ih{v) h"{u) 


AR4 = — [h{u — e) — h{u + £)]-j [h'{u — e) + h'{u + £)] = — 

where we substituted Eq. (E 3 ). Reverting to P2, we then obtain 

h{v) h'{u) 


{v — uY {v ~ 

8h'{u) 2h"'{u) 


v—u-\-s 


du 


AP4. 


f 

I2h'{u) 

h"'{u) ' 

L 

{v — uY 

{v — uY _ 


(ElO) 


dv. 


P 4 [h] = P J 


h"'iu) 


I-00 [{'>J ~ 'U-Y {v — uY 6(u — m) 2 J 
This leads to the main result of this subsection, which is that 


dv. 


l^P 

4 ! du'^ 


h{v) 


v — u 


dv = P 


1 


-00 (v - w)® L 


h{v) — h' {u){v — u) — h"'{u) 


{v — uY 

3 ! 


dv. 


(Ell) 


(E 12 ) 
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